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I. INTRODUCTION 


Direct Methanol Fuel Cells (DMFCs) are recently 
being attracted as an alternative power source to batteries 
for portable applications since they potentially provide 
better energy densities. However, there are two key 
constraints limiting the effectiveness of DMFC systems: 
crossover of methanol from anode to cathode and the 


sluggish kinetics of the electrochemical oxidation of 


methanol at the anode. 
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Abstract— A one dimensional (1-D), isothermal model for a direct methanol 
fuel cell (DMFC) is introduced and solved numerically by a simple finite 
difference scheme. By using numerical calculation, the model model can be 
extended to more complicated situation which can not be solved analytically. 
The model considers the kinetics of the multi-step methanol oxidation 
reaction at the anode. Diffusion and crossover of methanol are taken into 
account and the reduced potential of the cell due to the crossover is then 
estimated. The calculated results are compared to the experimental data 
from literature. This finite difference scheme can be rapidly solved with high 
accuracy and it is suitable for the extension of the model to more detail or to 


higher dimension. 


The crossover of methanol lessen the system 
efficiency and decreases cell potential due to corrosion at 
the cathode. The electrochemistry and transport processes 
in DMFCs are shown in Fig.l. Methanol is oxidized 
electrochemically at both the anode and cathode, however 
the corrosion current at the cathode does not create any 
useful work. A number of experimental and computational 


investigations have 


DMFCs [1-4]. 


reported methanol crossover in 
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co, CO: 
H-0 H 20 


CHOH 


o 
H20 i 


CH30H + H20 — CO2 + 6H* + 6em 


20, + 6H* + 6e > 3H,O 
CH:0H + H20 — CO2 + 6H* + Ge- 


50:2 + CH;OH — CO, + 6H* + 6e- 


Fig.1 Schematic illustration of a DMFC. 


There are several models have developed to predict 
the behaviour of direct methanol fuel cells, which is 
important in the design, operation and control. Among 
them, 1D model show the advantage of simple and fast 
calculation, which is suitable for real time simulation. 
Garcia et al. [5] presented a one dimensional, isothermal 
model of a DMFC to rapidly predict the polarization curve 
and goes insight into mass transfer happening inside the 
cell. The model was solved analytically. However, 
analytical methods have some drawbacks such as the 
limitation to some specific cases and difficulty to extend to 
more complicated situation. Therefore, in this current study, 
instead of using analytical method, the model is solved 
numerically using a simple finite difference scheme. 


One-dimension mathematical modeling of direct 
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methanol fuel cell 
The model which was developed in [5] is used in this study. 
The details are briefly discussed as follows. 
Assumptions. The model considers the 1D variation of 
methanol concentration across the fuel cell which includes 
anode backing layer (ABL), anode catalyst layer (ACL), 
and membrane. The schematic diagram of the layers 
considered in the model and several assumption illustration 
were presented in . The assumptions are detailed as 
follows 

1) Steady-state and isothermal operation. 

2) Variables are lumped along the flow direction 

3) Convection of methanol is neglected. 

4) Isothermal conditions. 

5) All physical properties, anodic and cathodic 
overpotentials are considered constant. 

6) Local equilibrium at interfaces between layers can 
be described by a partition function. 

7) All the reaction are considered as homogeneous 


reactions. 


membrane 


z 


0 55 65+ Sa Sat dat du 


Fig.2 Schematic diagram and concentration distribution of 


the DMFC layers 


The voltage of the cell is calculated as 
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in which, 
U o, and UMeOH are the thermodynamic 


equilibrium potential of oxygen reduction and methanol 
oxidation respectively 
Nc and na are the cathode and anode 
overpotentials, respectively 
Ou I Cell 
K 


represents the ohmic drop across the 


membrane. 
Anode backing layer - ABL (domain B) 
In this domain, the differential mass balance for methanol at 


steady state is 


B 
dN weon z 2 0 
dz 
(2) 
The methanol flux is the Fickian diffusion with an effective 
diffusivity Dg 
dc? 
B = MeOH ,z 
N MeOH,z —D,; 
dz 
(3) 


Combining Eq. (2) and Eq. (3), the distribution equation for 
methanol in ABL is 


2B 
d CMeOH,z 


dz” Ti 


(4) 
Anode Catalyst Layer - ACL (domain A) 
In this domain, there is a methanol oxidation reaction. 
Therefore, the differential mass balance for methanol at 


steady state is 


A 
dN MeOH,z _ MeOH 


dz M meon 
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(5) 


In which the molar rate of methanol consumption 


r, 
—MeOH _ is calculated from the volumetric current density 
M MeOH 
jas: 
meon _ J 
M moy OF 
(6) 


The current density is related to the concentration of 


methanol as ([6]) 


A 
: — qJ” kC weon a,n, FIRT 
J=a 0,ref A +A a,n, FIRT 
Cmeoy T A€ 
(7) 
f . f MeOH 
In which a is the specific surface area of the anode, T 0 se 


is the exchange current density, and k and å are constants 
The methanol flux is the Fickian diffusion with an effective 


diffusivity Da 


A 
NA -D dC eon.: 
MeOH,z TAa 
dz 


(8) 
Combining Eq. (5), Eq. (6) and Eq. (8), the distribution 
equation for methanol in ACL is 
2 A ; 
d Cmon: _ j 


D, — MOr = 
^ d 6F 


Membrane (domain M) 
The differential mass balance for methanol at steady state in 


the membrane is 


M 

AN weon z = 0 
dz 
Page | 302 


Nguyen et al. 


(10) 
The methanol flux in the membrane includes the diffusion 


and electro-osmotic drag as follows: 


M 
dC meon z 


dz 


I 
Cell 
+ Sueon F 


(11) 
In which Dm and meon are the effective diffusion in 
membrane and the electro-osmotic drag coefficients of 
methanol, respectively. 
Combining Eq. (10) and Eq. (11), the distribution equation 


for methanol in membrane is 


2 M 
d CMeOH,z 


Pu dz? 


=0 
(12) 

Boundary condition: 

At z=0 (the interface between the flow-channel and anode 

backing layer), there is no mass resistance. Therefore, the 

concentration is given by the bulk concentration of the flow 


as: 


B = 
C MeOH ,z=0 T Cbulk 


(13) 
At z= op (the interface between ABL and ACL), there are 
two conditions. First, the local equilirium of the 
concentrations between two domains is given by a partition 


coefficient Kı as 


A _ B 
CMeOH .2=5, ~ K ‚Cheon 2-8, 


(14) 
Second condition is the equality of fluxes between two 


domains (ABL and ACL) 


B a. 
N MeOH ,2=5, 7 N MeoH z=5, 


(15) 
At z= ôs + da (the interface between ACL and membrane), 


there are two conditions. First, the local equilirium of the 
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concentrations between two domains is given by a partition 


coefficient Ky as 


M =K g^ 
C MeOH ,z=6,+8, T +$ HC MeOH ,z=6,+6, 


(16) 
Second condition is the equality of fluxes between two 


domains (ACL and membrane) as 


A _ ag 
N MeoH z=5,+5, _ N MeOH 25,465, 


(17) 
At z= 6g+ ĝa + 6m: All the methanol crossing the membrane 
is assumed to consume immediately at the cathode, result in 
a zero concentration at the membrane/ cathode-layer 


interface. Thus, 


a =0 
C MeOH ,z=6,+6,+6, B 


(18) 
Finite difference scheme and overpotential calculation 
The spatial independent variable z in the three segments (0, 
Op), (dp, 63+ ôa), (Opt ĝa, Opt ôa + ôm) can be discretized 


into np, na, Nu Subdivisions, respectively, as 
O= <2 <.< Z, =ð; 
(19) 
p= <z << Ze =f t+ Oy 
(20) 
Op tO, = BY Sa <..< 2p = 63 +6, 4 Oy 


(21) 
In each segment, note that the length of subsegment is 
equal to Azp, Aza, Az, respectively. 
Governing equations 
Inside the domains (ABL, ACL and membrane), the 
second derivatives in the governing equations are 
discretized using central difference formulae. The details 
are as follows 


In ABL region, equation (4) is discretized as: 
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B B B 
CMeOH z+Az ~ ZC MeOH ,z + CMeOH 2-A: =, 


(Ace) 


B 


(22) 


B B B _ 
Cueou ist ~ 2CMeoH i + Cmeon i = 9 


(23) 


In ACL region, equation (4) is discretized as: 


A A ich 
C MeOH ,z+Az — ©MeOH,z T ©MeOH,z+Az 


A 2 
(Az, ) 
A 
q[ Mou KC Meon z e2 FIRT 
O,ref A +A a,n, FIRT 
CmeoH,z T ^E 
6F 
(24) 


A 2 
(Az, ) 
A 

pMeoH Keon i e2 FIRT 
a O,ref A +A a,n, FIRT 
Cueou i T ^e 
6F 
(25) 


In membrane region, equation (12) is discretized as : 


M M M 
C MeOH ,z+Az ~ ZCMeOH,z + CMeOH zz -0 


(Azu y 


M 


(26) 


M M M _ 
Cueou ist ~ 2CMeoH i + Cmeon i = 9 


(27) 
Boundary conditions 
The first derivatives in boundary conditions are 
approximated using forward difference formulae as 
follows: 


At the left interface, using the forward scheme: 
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AC yeon,c _ Cmeou,c+hc ~ CMeon z 


dz Az 


(28) 


At the right interface, using the backward scheme: 


dC yon z _ ĈMeOH,z ~ Ĉ MeOH ,z-Az 


dz Az 


(29) 
Concentration profile 
After discretization, a system of equations for the 
concentration of methanol is obtained. The system is 
solved using simple iteration method to find the 
concentration profile of methanol. 
Anode overpotential 
From the concentration profile, the cell current can be 


estimated as: 


6, +6, 


A 
Pe cede ql MoH kC meon 
cell — 0,ref 
5, 


a,n, FIRT 


a,n, FIRT 


A 
Cmeoy + 2€ 


(30) 

In which ya is assumed to be constant. The integration is 
numerically calculated using trapezoidal rule. Because 7 
is also included in calculation of concentration profile, an 
iteration is required to find appropriate ya for a given value 
of Icen. 

Cathode overpotential 

Tafel kinetics with first-order oxygen concentration 
dependence is used to estimate the oxygen reduction at the 


cathode. 


I 


cell 


+ Irak = 0, ref 


(31) 
In which /icak is the leakage current density due to the 
oxidation of methanol crossing the membrane. The leakage 


current density can be estimated as 
_ M 
Deak = 6FN weon, 


(32) 
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In which NG cous is estimated from Eq. (11). Then, Eq. 


(32) is used to obtain yc for a given value of Icen. 


Table /. 
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After the anode and cathode overpotentials are known, the 
Vceu for a given value of Icen is calculated using Eq. (1). 


The parameters used in the model are summarized in 


Table 1 Model parameters 


Parameter Value 
a 1000 cm"! 
Da 2.8 x 10°exp(2436(1/353-1/T)) cm?/s 
Dg 8.7 x 10% cm?/s 
Dm 4.9 x 10°°exp(2436(1/333-1/T)) cm?/s 
Do 9.425 x 103exp(33570/R(1/333-1/T)) A/cm? 
li 4.222 x 103exp(73200/R(1/333-1/T)) A/cm? 
Kr 0.8 
Ku 0.8 
k 7.5 x 10+ 
T 343.15 K 
Umeou 0.03 V 
Uo2 1.24V 
Oa 0.52 
Oc 1.55 
Oa 0.0023 cm 
ôB 0.015 cm 
ÔM 0.018 cm 
K 0.036 s/cm 
A 2.8 x 10° mol/cm? 
¢MeOH 2.5XMeOH 
Il. RESULTS AND DISCUSSIONS the end of the curve is quite high. The disagreement could 


The simulation results of the polarization curve for DMFC 
at different concentrations of the bulk flow are shown in 
Fig.3. The calculation results well agree with the 


experimental data report in [5]. However, the difference at 
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be due to the assumption that the methanol electro-osmotic 
drag coefficient is a constant value. It is better to calculate 
the electro-osmotic drag coefficient at each point, especially 


at the end of the curve. 
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0.8 


0.4 


cell 


Fig.3 Model predictions for different methanol concentrations 


Fig.4 shows concentration profiles across the anode and membrane obtained by the model for the four concentrations at 
15 mA/cm’. 


0.5 


o o > 
io w P 


Methanol concentration (mol/L) 
> 


0.0 = 
0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 


z (cm) 


Fig.4 Concentrations profiles for different methanol bulk concentrations 


MI. CONCLUSIONS parameters from literature, the calculation results well 
In this study, a finite difference scheme were sucessfully agree with experimental polarization curve. The scheme 


applied to solve the one-dimensional, isothermal model of also is applicable in the estimation of concentration 


a DMFC. Using reasonable transport and kinetic profiles in the anode and membrane as well as predicting 
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the methanol crossover. The computation time is fast 


enough for real time application. 
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